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We have developed a new numerical scheme to solve r-mode oscillations of rapidly rotating poly- 
tropic stars in Newtonian gravity. In this scheme, Euler perturbations of the density, three compo- 
nents of the velocity are treated as four unknown quantities together with the oscillation frequency. 
For the basic equations of oscillations, the compatibility equations are used instead of the linearized 
04 , equations of motion. 

j^Q' By using this scheme, we have solved the classical r-mode oscillations of rotational equilibrium 

' sequences of polytropes with the polytropic indices TV = 0.5, 1.0 and 1.5 for m = 2, 3 and 4 modes. 

Here m is the rank of the spherical harmonics Y t m . These results have been applied to investigate 
evolution of uniformly rotating hot young neutron stars by considering the effect of gravitational 
radiation and viscosity. We have found that the maximum angular velocities of neutron stars are 
around 10-20% of the Keplerian angular velocity irrespective of the softness of matter. This confirms 
the results obtained from the analysis of r-modes with the slow rotation approximation employed 
by many authors. 
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I. INTRODUCTION 
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About 110 years ago, Bryan analyzed oscillations of self-gravitating and spheroidal liquids and obtained the period 
equations for oscillation modes whose frequencies are proportional to the angular velocity Q). These modes are 
' similar to those which were investigated for rotating realistic compressible stars and named r-modes by Papaloizou & 
Pringle | (see also e.g. |fl). 

It is only very recently that the r-mode has been rediscovered from a standpoint of its strong instability and 
intensively studied. For sufficiently rapidly rotating stars whose oscillations are coupled to gravitational radiation, 
instability induced by gravitational radiation was discovered by Chandrasekh ar pj and by Friedman & Schutz || . 
• i-H Thus it is called the CFS instability (for reviews, see e.g. §-0). Ander sson |1 lfl has pointed out that the r-mode 
oscillations suffer from the CFS instability and play a very important role in a first few years after formation of 
\ neutron stars. 

In subsequent studies by many authors (e.g. |12-17|), it has become clear that the r-mode oscillations make newly 
born neutron stars significantly unstable by gravitational radiation emission for a temperature range around 10 9 K. 
For the CFS instability of the f-modc oscillations, the effect of viscosity stabilizes neutron stars almost all range of 
the temperature. Consequently the angular velocity would be affected only slightly and would become 90-95% of the 
Keplerian angular velocity. 

Contrast to the CFS instability of f-modes, the CFS instability of r-modes will not be stabilized for a certain range 
of the temperature T of neutron stars, i.e., roughly around T ~ iq8~io K, even if there exists viscosity. This may be 
explained by the weakness of bulk viscosity damping because the volume expansion rate of the r-mode oscillation is 
small. Also the critical rotational frequency of the star from which the instability sets in is zero, in contrast to the 
case of the f-modes. It leads to the enhancement of the instability as well. Due to the effect of the r-mode instability, 
neutron stars which have been born with T ~ 10 11 K will lose most of their angular momentum in only one year 
and settle down to slowly rotating configurations. This time scale is determined from that of cooling of newly born 
neutron stars down to T ~ 10 9 K. Thus there exist critical angular velocities for neutron stars above which no stable 
configurations exist for each value of the temperature. 

This sets a severe limit to the scenario of evolution of newly born neutron stars, in particular, to that of millisecond 
pulsars. From the observational point of view, the rotational periods of neutron stars were considered to be rather 
long (see e.g. |lq,n^]), although from the angular momentum considerations of the progenitors of neutron stars they 
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could be shorter than ~ 10 ms. Therefore, some mechanism which sets limit to the rotational periods of neutron stars 
has been required and looked for. From the significant feature of the r-mode instability, it could be the promising one 
which explains the observational results. 

However, we should note that almost all results mentioned above have been obtained from stability analysis of slowly 
rotating stars except two works by Bryan 111 and Lindblom & Ipscr |l4j. For example, Lockitch & Friedman |l5| ] 
have introduced the slow rotation approximation to get eigenfrequencies and eigenfunctions of the classical r-modes 
and generalized r-modes for polytropic stars (about generalized r-modes, see also |l^]), and Andersson, Kokkotas, 
& Schutz [fl6| have also analyzed slowly rotating equilibrium configurations to estimate the spin evolution of young 
neutron stars. As for the results of Q| and those of |l4|], Maclaurin spheroids were analyzed. Since realistic stars are 
not incompressible fluid and newly born neutron stars might be rotating with very short periods, it would be desirable 
to investigate r-mode oscillations of rapidly rotating and compressible stars. 

In the previous Letter paper J2(]] , we have shown some representative results of the classical r-mode oscillations of 
rapidly rotating polytropic stars with the polytropic index N = 1. Those results have been obtained by developing a 
new scheme to handle r-mode oscillations of rapidly rotating polytropic stars in Newtonian gravity. In this paper, we 
will explain the new scheme in detail and show classical r-mode oscillations of rotational equilibrium configurations 
for a wide range of N and m. 

The numerical scheme is an extended version of the one which was used to analyze f-mode instability by Yoshida 
& Eriguchi |21|. After computing r-mode oscillations for equilibrium sequences of several rotating polytropes, we 
will calculate time scales of evolution due to dissipation by taking account of the effect of viscosity and gravitational 
radiation and obtain critical curves for the r-mode instability in the temperature - rotational frequency plane by 
following the analysis method proposed by Ipser & Lindblom p2). In the same plane, we can also draw "evolutionary 
curves" due to cooling by using the standard modified URCA process [ p3| or other mechanisms. 

Our newly developed scheme for the oscillations of rotating stars will be explained in Section 2. In Section 3, we will 
show our computational results for the classical r-mode oscillations of rotating polytropes. Time scales of evolution 
due to gravitational radiation as well as viscosity will be calculated by using the eigenfrequencies and eigenfunctions 
of r-mode oscillations. By using these time scales we will be able to get critical curves for the r-mode instability 
and the 'evolutionary' curves for hot young neutron stars. Finally, in Section 4, we will discuss the uncertain issues 
included in our analysis, and some other important problems which we will have to treat carefully in the future. 



II. SCHEME OF LINEAR PERTURBATIONS FOR UNIFORMLY ROTATING POLYTROPES 



A. Equilibrium states 



In this paper we will treat uniformly rotating polytropes in the framework of Newtonian gravity. Axisymmetric 
equilibrium configurations of rapidly rotating Newtonian stars are obtained by solving the hydrostatic equation, 
Poisson's equation and the equation of state as follows: 

— Vp = - V0o + r sin 9n 2 e R , (1) 
Pa 



A^o = 4ttGpo , (2) 

po = Kpl + * , (3) 

where po, po, </>o, Q, en, G, K and N are the density, the pressure, the gravitational potential, the angular velocity, a 
unit vector in the R direction of the cylindrical polar coordinate (R, (p, z), the gravitational constant, the polytropic 
constant, and the polytropic index, respectively. Through this paper, the spherical polar coordinates (r, 9, (p) are used. 

In order to handle properly the boundary condition of the potential at infinity, the potential is transformed into the 
integral representation by using the Green's function. These equations are solved by applying the method developed 

In actual computations we have obtained three polytropic equilibrium sequences with N = 0.5, 1.0, and 1.5. The 
mesh number in this paper is (r X 6) = (32 x 11). Here equidistant mesh sizes in r- and (^-directions are used. 
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B. Linear perturbations of uniformly rotating and axisymmetric polytropes 



Once equilibrium configurations are obtained, we can proceed to analysis of linear perturbations of equilibrium 
states. Basic equations for the perturbed states can be written as follows. We will express Eulcrian perturbations by 
using 5, say Sg for the physical quantity g. Since in this paper rotating polytropes will be employed, the linearized 
equations need to be treated separately for N — and N ^ cases. 



1. Basic equations for N / polytropes 

For TV 7^ polytropes, the linearized equations of the continuity and the linearized equations of motion are written 
as follows: 

1 dSpi dSui 2 . 1 dpo r 1 dSvi 1 1 dp a . 

+ -rr 1 + Sui + ^-5 Ul + — - + kf-Sv! 



po dt dr r po dr r 89 r po 80 

1 „<- ^ 1 85 pi 1 88wi 
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Here perturbed quantities of the density three components of the velocity, (u, v,w) in the spherical coordinates, the 
gravitational potential are expressed as Spi,5ui,5vi,Swi and 5<j>i, respectively. We have made use of the relation 
between the perturbed pressure, Spi, and the perturbed density as follows: 

5p 1 =(l + ^K p y N 5p 1 . (8) 

Since the unperturbed states are axisymmetric and stationary, all perturbed quantities can be expressed by using 
the expansion as follows: 



Spi{r,0,ip,t) 


= ^ exp(i(crt — 
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5<j)i(r,9,ip,t) 


= ^ cxp(i(crt — 

m 


mtp))4> m (r,9) , 



(9) 



where cr, p TO , u m , u m , w m and m are the oscillation frequency and the expansion coefficients of corresponding quanti- 
ties, respectively. 

When we define the physical quantities 
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Sp = Prn , 

Svg = in 



5w v = w 



m 5 
m ; 



and apply the above expansion, the perturbed equations become as follows: 
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(11) 
(12) 

(13) 

(14) 



2. Basic equations for N — polytropes 



For N — polytropes, we can follow the same procedure as that for N ^ polytropes by choosing the pressure 
perturbation as one of the unknown quantities instead of the density perturbation. The final equations can be written 
as follows: 
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(15) 
(16) 
(17) 
(18) 



3. Compatibility conditions 



Instead of solving the basic equations derived above, we make use of the compatibility relations among the equations 
of motion, in particular, the compatibility equations between Eqs. (|l^) and ( |l4| ) and between Eqs. (|lj) and ( |l3| ) for 
N 7^ models or those between Eqs. (|lf) and (|l|) and between Eqs. © and (0) for N = models. They can be 
written both for N ^ and N = cases as follows: 
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4- Perturbed gravitational potential 



The perturbed gravitational potential can be expressed by using the integral form as follows: 

o/„:-a> {n-m)\ 



= -4ttG^ / dO'sm 



(n + to)! 
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(21) 



where r s (9) is the surface radius of the equilibrium configuration and Sr s is the change of the surface. Here, functions 
/„(r, r') are defined as 



1 ir 1 



(22) 



and P™'s are associated Legendre polynomials. 

It is noted that the boundary conditions for the gravitational potential, i.e. regularity condition throughout the 
space and the flatness at infinity, are automatically included in this integral representation. 



5. Boundary conditions on the stellar surface 



Since the perturbed gas must flow along the perturbed surface, the boundary condition for the flow velocity on the 
surface (for N ^ configurations) or that for the surface (for N — configurations) can be expressed as follows: 



for N 7^ polytropes, or 



^-Su r + - ^-5v e + (mil - a)Sp = 
or r oO 



Svg dr B 

5u r H — - + (mil — a)dr s = 

r s dO 



(23) 



(24) 



for N = polytropes. 

For N = polytropes, since the unknown quantity 5r s has to be solved in addition to other physical quantities, we 
need to include the following relation which is satisfied on the surface: 



Sp s 



( ^El 
{ Or 



Sr s = po 



- r s sin 2 9il 2 



5r s . 



(25) 



6. Solving scheme 



Our basic equations and boundary conditions are Eqs. (11), (|l4|), (19), (p0|), and ( |23] ) for N ^ polytropes and 
Eqs. (HI), (H|), ©, ©, ©, and (||) for N = polytropes, although we will not show the results for TV = 
polytropes in this paper. 

In order to handle the surface boundary conditions as precisely as possible, we introduce surface-fitted coordinates 
as follows (see e.g. p3-p6|): 



r s (0) 



(26) 
(27) 
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We choose equidistantly spaced mesh points on the r*- and 6**-coordinates. The basic equations expressed in the 
new coordinate system are discretized on these mesh points. The derivatives with respect to r* and 9* are replaced 
by the central difference scheme with the second order accuracy for the most mesh points except for those near the 
surface and the equatorial plane. The derivatives with respect to r* near the surface and those with respect to 9* 
near the equator are replaced by difference schemes with higher order accuracy. 

Instead of solving the eigenvalue problem directly, we add one more condition to the basic equations of the system 
and solve for all the unknown quantities simultaneously by using the Newton-Raphson iteration scheme. Since we 
focus only on the classical r-mode in this paper, we adopt the following condition as the additional equation: 



We have computed equilibrium sequences of uniformly rotating polytropes with N = 0.5, 1.0 and 1.5 by starting 
from spherical models to terminal models which are configurations just before shedding mass from the equatorial 
surface. By using these equilibrium models, classical r-mode oscillations for m = 2, 3 and 4 modes have been solved 
and eigenfrequencies and eigenfunctions have been obtained. In the actual numerical computations, the mesh number 
in (r*, 9*) coordinates is the same as that used in obtaining the equilibrium configurations, i.e. (r* x 9*) = (32 x 11). 



In Figs. 1-3, the eigenvalue of the m = 2 r-mode is plotted against the dimensionless angular velocity of the 
equilibrium configurations for N = 1.0,1.5. and 0.5 polytropes, respectively. In these figures, the eigenfrequency a 
is normalized by using the angular velocity f2 and the angular velocity is normalized as f2/ \/ATrGp c . Here p c is the 
central density of the star. 

Since classical r-mode oscillations of polytropes have been investigated by using the slow rotation approximation 
to the third order of the angular velocity in |17|| , we can compare our results with theirs. As seen from these figures 
except Fig. 3, the relative errors between the numerical results and those of the approximate results are 0.4% for 
the slowest configurations. It is natural that for rather rapid rotational models difference becomes larger because the 
assumption of slow rotation is violated. 

Fig. 4 shows the eigenfrequencies of N = 1.0 rotating polytropes for different values of m. From this figure, it is 
clear that the behaviors of the eigenfrequency along the rotational sequence for different values of m are almost the 
same. This is also found for N = incompressible configurations, i.e. for Maclaurin spheroids 0]. In Fig. 5, the 
eigenfrequencies of N — 0.5, 1.0 and 1.5 polytropes for m = 3 perturbations are plotted. 

In Figs. 6-7, we have plotted the profile of the perturbed quantities, (Su r , Svg, Sw^, Sp), against the radial distance 
normalized by the stellar radius, i.e. r*. Figs. 6 and 7 display the distributions of the r-mode oscillation of the 
polytropes with N = 1.0 for r p /r e — 0.99 and r p /r e = 0.71 models, respectively. Here r p and r e are the polar radius 
and the equatorial radius, respectively. 



In order to analyze the stability due to dissipative mechanisms such as gravitational radiation and viscosity, there 
are several schemes to include the dissipative forces. A direct method is to solve the perturbed equations of motion in 
which the forces due to gravitational radiation and viscosity are taken into account [|2lj . Another method is to evaluate 
the time scales of the system change due to gravitational radiation and viscosity |22j | (see also [p|-|To|). Although two 
methods are essentially the same, the formulation of the latter method is very simple. Thus in this paper we will 
follow the Ipscr-Lindblom method. 

In their method, the time derivative of the energy of the oscillation mode plays an important role. The energy of 
the mode, E(t), as measured in the rotating frame can be expressed as 



Svg — 1.0 , on the equatorial surface . 



(28) 



III. NUMERICAL RESULTS 



A. Eigenfrequencies and eigenfunctions of the r-mode oscillations 



B. CFS instability 




(29) 



G 



where * represents complex conjugation. From this definition and by using the perturbed continuity, the perturbed 
equations of motion with the dissipative forces and the perturbed Poisson's equation, we can get the following 
expression for the time derivative of the mode energy in the rotating frame: 



dE _ dE gr dE s dE b 
dt dt dt dt 



(30) 



where 

dEg r 

dt 

dE s 
dt 

dEt 
dt 



-a{a-mn) ^ N t a 21 (\8D lm \ 2 + |<5J; m | 2 ) , (31) 

l>2,l>m 

- f 2 V 6a ab Sa* ab d 3 x , (32) 



(\S@\ 2 d 3 x. (33) 



Here Ni, D[ m , Ji m ,r], £, a ab and are a coupling constant, the mass multipole, the mass current multipole, the 
coefficient of shear viscosity, the coefficient of bulk viscosity, the shear and the expansion, respectively. They are 
defined as: 

= 4ttG (Z + !)(/ + 2) 

c 2/ +U(Z-l)[(2l + l)!!] 2 ' 1 ' 

D lm = Jr l pYr i d 3 x , (35) 

J lm = H J_ f r l (pv) ■ (fx VY t * m )d 3 x , (36) 
cl + 1 J 

°ab = V (a u b) - — S ab , (37) 

6 ee V c w c , (38) 

where c is the speed of light and a round bracket in the subscript means symmetrization among indices. 

For slowly rotating Newtonian stars, the eigenfrequency to the first order of the angular velocity is expressed as 
II 

fJslow = y 1 ~ TQ + rj ) ' (39) 

From this eigenfrequency for slowly rotating Newtonian stars, the coefficient of the energy dissipation rate due to 
gravitational radiation Eq. (p^), i.e. the first term in the righthand side of Eq. (pOh, becomes positive as follows: 



a slow (a slow - mO) = ^-^tt^^^ > . (40) 

r(t + \y 



Consequently, since the time derivative of the energy dissipation due to gravitational radiation becomes positive, 
gravitational radiation makes the system unstable by growing the perturbation. 

On the other hand, as seen from the second and third terms in the righthand side of Eq. (|30"|), viscosity works as 
a stabilizing factor for the system. Therefore, if we compare the time scale of gravitational radiation and that of 
viscosity, we will be able to conclude the stability of the system. 

The time scale of the system change due to dissipation can be defined as the inverse of the imaginary part of the 
mode eigenfrequency, r r as follows: 

l_ = E gr + E s +E b ^ J_ + j_ + 1_ ^ 

where ' = dj dt. Positive values of r" 1 imply that the effect of gravitational radiation emission overcomes the stabi- 
lization effect due to viscosity. In other words, the system is in an unstable state. Therefore, the critical condition 
which divides stable and unstable configurations can be defined by the following equation: 



7 



(42) 



In order to evaluate the time scales we need to know the viscosity coefficients at high density regions. As for the 
origin of the viscosity, we adopt the same mechanisms as p2fl . The shear viscosity is owing to neutron-neutron or 
electron-electron collisions in the neutron matter when the temperature of the star is not so high. On the other hand, 
for high temperature regions, the bulk viscosity which is due to the energy generation caused by compression of the 
neutron star matter because the time scale to balance between /3-decay and inverse /3-decay is longer than that of 
the dynamical perturbed motion of the star is dominant. The coefficients of each viscosity are written as follows 

Vnn = 2.0 x 10 18 p% 4 T~ 2 g cm-V 1 , (43) 

77ee = 6.0 x 10 18 p 2 15 Tg 2 g cm-V 1 , (44) 

C = 6.0 x 10 25 p? 5 (o- - mQ)~ 2 T$ g cm^s" 1 , (45) 

(46) 

where 

P 15 = -mis P =3 ' ( 47 ) 

10 ia g cm 6 

For lower temperature regions, we need to consider only electron-electron collisions for the coefficient of shear. 



C. Evolution of hot young neutron stars 



In this paper we will choose the polytropic constant K so that the mass and the radius of the star become 
M = IAMq and R w 12.5km in the spherical limit. We will call this star a canonical neutron star. If we specify the 
polytropic index, the value of the polytropic constant K is uniquely determined. 



1. Critical curve for the r-mode instability 

The critical condition (^) depends on the structure of the unperturbed configuration, the eigenfunctions of the 
perturbed state, the eigenfrequency and the temperature. Since we assume the polytropic relation for the canonical 
neutron stars and consider only classical r-modes, the structure of the equilibrium state, the eigenfunctions and the 
eigenfrequency are all determined by specifying the angular velocity in addition to the polytropic index N and the 
mode number m. Thus the critical condition can be considered to be a function of the angular velocity and the 
temperature as well as N and m. In other words, once we specify the values of N and m, the critical condition can 
be expressed as a critical curve on the temperature - rotational frequency plane. 



2. 'Evolutionary' curve for the stellar spin 

It is not easy to follow realistic evolution of neutron stars by considering the effect of gravitational loss and viscosity 
as well as the effect of cooling due to neutrino loss. Thus we will choose a different approach. 

Once cooling mechanisms are specified, we can estimate cooling time scales of the neutron star. In this paper, we 
assume that the standard modified URCA process is dominant in the initial stage of the neutron star formation [E3] . 
For this process, the effective cooling time scale is a function of the stellar temperature and defined by 



dhiTg 
dt 



£. (49) 



where t c is a cooling time for the modified URCA process, typically ~ ly. 

Since the neutron stars evolve due to radiation of gravitational wave and neutrino loss, we can define a critical state 
where the time scales of two mechanisms are the same as follows: 
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— ''"cool 



(50) 



Since this state is also a function of the rotational frequency and the temperature if N and m are specified, we can 
draw an 'evolutionary' curve of the stellar spin, whose meaning is made clear in the following discussion, on the 
rotational frequency-temperature plane. 

3. Spin evolution of canonical neutron stars 

In Figs. 8-12, critical curves for the classical r-modes instability are drawn in the rotational frequency, /, and the 
temperature plane for different values of N and m. Above these solid curves, the instability due to gravitational wave 
emission overcomes the stabilizing effect of viscosity and consequently neutron stars become unstable. 

In these figures, the evolutionary curves of the stellar spin are also shown by dashed curves. Above these curves, 
cooling time scale of each cooling mechanism is longer than the instability time scale. It implies that in this region 
the neutron stars evolve nearly downwards by losing the angular momentum but keeping the temperature constant. 
In the region below these curves, the stars sweep leftwards due to neutrino loss by keeping the rotational frequency 
almost constant. 

Suppose we fix a cooling mechanism of the star. It is certain that neutron stars are born with high temperature 
with T ~ 10 11 K. Concerning the rotational frequency of the newly born neutron stars, they may not rotate so rapidly 
compared with the rotational frequency of the Kepler limit [fl8| , [l9[ but may rotate with certain rotational frequencies 
above the minimum value of the critical curve. It implies that such newly born neutron stars are clearly located below 
the critical curve at the right side of the T—f plane. Since this region is also below the evolutionary curve, cooling 
stars move almost horizontally to leftwards on this plane. 

During the evolution due to cooling, the stars come to and cross the critical curve. Once they come into the regions 
above the critical curve, the r-mode instability due to gravitational radiation loss starts working significantly. As 
a result, the stars would evolve downwards vertically. However, the stars are still located below the evolutionary 
curve and the time scale of cooling is smaller than that of the instability. Therefore, due to this cooling effect, the 
stars move to leftwards almost horizontally on the plane. Eventually the stars go across the evolutionary curve. As 
soon as they cross the curve, the instability works efficiently and the stars move downwards vertically. If they cross 
the evolutionary curve once more, the cooling effect becomes dominant again and the stars will evolve horizontally. 
Consequently we can consider that the stars evolve almost along the evolutionary curve of the cooling mechanism 
being considered. In this sense we can regard the evolutionary curves as curves along which the spin evolution of a 
cooling star takes place. 

As seen from Figs. 8-12, the evolutionary curves are located above the critical curves. This means that the evolution 
mentioned above will come to an end because they go out of the unstable region of the gravitational radiation and 
reach the lower left part of the T—f plane. Consequently, the neutron stars will settle down to the states with the 
rotational frequencies corresponding to the minimum values of the corresponding evolutionary curve for each cooling 
channel. 

In Figs. 8-10, we show the critical and evolutionary curves for m = 2, m = 3 and m = 4 modes, respectively. From 
these figures, the minimum values of critical curves are about 100Hz, 200Hz, and 300Hz for the m = 2, m = 3, and 
m = 4, respectively. Therefore we can see that the instability of the m = 2 mode dominates other modes with higher 
values of m. 

In Figs. 11-12, we show the critical and evolutionary curves of polytropes with N — 0.5 and N — 1.5 for the m = 2 
r-mode. These results are almost the same as the result of N = 1.0 polytropes. It means that the classical r-mode 
instability depends little on the compressibility. This agrees with the result of jl3| in which they conclude that the 
instability is insensitive to the equation of state. In particular, the smallest values of the rotational frequency of the 
critical curves seem to be determined uniquely irrespective of the polytropic index as far as the value of m is the same. 
Thus we can conclude that neutron stars will settle down to states with a unique rotational frequency, i.e. around 
100Hz. 

As a summary, neutron stars which are born with high temperature (lO^K) and rather rapid rotational frequency 
(a few times of 100Hz) will be spun down to about 100Hz by the effect of the gravitational instability on a quite short 
time scale. 



IV. DISCUSSION 

In this paper we have developed a new numerical method for solving the r-mode oscillations of rapidly rotating 
compressible stars. By using the obtained eigenfrequencies and eigenfunctions, we have calculated the condition for 
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occurrence of the gravitational radiation induced instability. By combining the cooling of the neutron stars and the 
instability, we have confirmed the evolutionary scenario of hot young neutron stars proposed from the slow rotation 
approximation, in particular, considerable spin down of the neutron stars with rather high rotation rates of initial 
spins on a very short time scale (e.g. [jl3|j29| ). The shortness of the spin down time scale implies that this mechanism 
would, effectively, result in "birth of slowly rotating neutron stars" because it would be difficult to observe the very 
initial stage with rather high rotational frequencies during one or so years just after the formation of neutron stars. 
This might be the reason of observational data of nonexistence of pulsars with very short initial spin periods Jl8],|l9| . 

It should be noted that our analysis has been carried out under some restricted assumptions. First, we have 
formulated the problem in the framework of Newtonian gravity. In order to get quantitatively correct values we 
have to use general relativity, although it is considerably difficult to solve oscillations of not only the matter but also 
the gravitational field. At the moment, no one has succeeded in treating the perturbations of the gravitational field 
of axisymmetric configurations except obtaining neutral points for f-mode oscillations p2| . Concerning the r-mode 
oscillations in general relativity, there arises another problem. Since there is the frame dragging, the effective angular 
velocity cannot be uniform any more. Some authors suggest that this may make the problem a singular eigenvalue 
problem || (see also @]). 

Second, we have to treat the realistic equation of state for the neutron star matter instead of simplified equations 
of state such as polytropes. However, as has been known from the slow rotation approximation J^-Q] and as we have 
shown for rapidly rotating models in this paper, there arise little differences of the oscillations of polytropes with 
different polytropic indices. Therefore, even if we treat the realistic equation of state, the situation would not change 
drastically as far as the classical r-modes are concerned. 

Third, there is some uncertainty about the cooling process. Although the modified URCA process has been taken 
as the standard mechanism of the neutron star cooling, there are some other mechanisms which would contribute to 
cooling of neutron stars. Therefore, we have computed evolutionary curves by introducing several possible cooling 
mechanisms. Fig. 13 shows such curves for several cooling processes: the neutrino bremsstrahlung, the quark (3 decay, 
the pion condensation, and the direct URCA cooling model, beside the modified URCA process. In our calculations, 
the time scale of the direct URCA process is taken from |}(]] and those of other processes are taken from HU (see also 
plf). As seen from this figure, other cooling mechanisms give obviously different evolutionary curves or "evolutionary 
paths". Concerning the minimum values of the evolutionary curves, however, differences among different cooling 
processes are rather small. In other words, the final rotational frequencies depend weakly on the cooling mechanisms. 
Consequently the final rotational frequencies would be around 50 ~ 180Hz as far as the scenario mentioned above 
works for the hot young neutron stars. 

Fourth, the assumption that neutron stars rotate uniformly is also uncertain. Since we are dealing with the very 
beginning of the neutron star formation stages, newly born neutron stars may rotate differentially, although they will 
soon settle down to uniform rotation on a time scale of several years. The r-mode oscillations of differentially rotating 
stars are another problem. There may arise a corotation point within the star and the similar problem as that in 
general relativity needs to be solved. 

In addition to these problems, we have to discuss the following crucial issues which may play a more important role 
in realistic formation and evolution of neutron stars than the above mentioned things. 

The superfluidity is one of the most important but poorly understood physics about the realistic neutron stars. Some 
authors (e.g. p5| ) have suggested that the effect of the superfluid might suppress the r-mode instability completely in 
the low temperature regions, while other authors have assumed that this effect might not be strong enough (if]] (see 
also ]l6|). In this paper, we assume that the superfluidity will not work effectively. In fact, since the realistic neutron 
star matter cannot be examined by any present experiments on the earth, the only way to estimate this effect may 
be the comparison between the theoretical works and the observational data. 

The existence of a solid crust may be another important factor for the r-mode instability. Bildsten & Ushomirsky |36| ] 
have suggested that the existence of the solid crust whose thickness is ~ 1 km would work as a stabilizing factor of 
the r-mode instability and that the minimum rotational frequency might be 40% of that of the Kepler limit. However, 
it is not clear whether or not the solid crust has formed only in a year or so from the formation of the neutron star. 
Furthermore, it is uncertain about the transition region from the fluid to the solid. Thus we have not taken the effect 
of the crust into our considerations in this paper. 

The most important factor concerning the neutron star evolution may be the presence of the magnetic field. There 
is a high possibility that hot young neutron stars might spin down by the effect of the magnetic field and/or that 
neutron stars might be born with rather slow rotation rates because the angular momentum is transferred by the 
effect of the magnetic field from the core to the envelope during the evolution of the massive stars (e.g. [ 
Inclusion of magnetic fields may change the r-mode characteristics itself. For instance, magnetic braking may enhance 
the instability by increasing the amplitude of the mode |38| . Also suggested is that the CFS instability by the Alfven 
wave radiation may exist for stars with strong magnetic fields. On the other hand it is argued |39| ] that the drifting 
fluid motion produced by the r-mode may generate strong toroidal magnetic fields. This may suppress the growth 
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of the r-mode amplitude. Therefore, the realistic evolution of the neutron stars have to be studied by including the 
r-mode instability as well as the effect of magnetic fields in the future. 
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FIG. 1. The eigenfrequency normalized by the angular velocity, a/Q, of the m — 2 r-mode is plotted against the dimen- 
sionless angular velocity, Q,/ \/ AnGp c , for N = 1.0 polytropes. Solid and dashed curves show our numerical result and the 
approximate result obtained in jnj], respectively. The relative difference between them, at the slowest rotational frequency, i.e. 
0,/yJA-KGpc w 0.02, is 0.4%. 
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FIG. 2. Same as Fig. 1 but for TV = 1.5 polytropes. The sequence terminates around fl/ ^4nCp c ~ 0.15 because it is difficult 
to obtain smooth eigenf unctions for more rapidly rotating configurations. 

FIG. 3. Same as Fig. 1 but for N = 0.5 polytropes. Since there are no results of the slow rotation approximation in Jl7| ], 
only our results are shown. 

FIG. 4. Same as Fig. 1 but for different values of m, i.e. m = 2 (solid curve), m — 3 (dashed curve) and m = 4 (dotted 
curve) . 

FIG. 5. Same as Fig. 1 but for m = 3 and for different values of TV, i.e. TV = 0.5 (dashed curve), TV = 1.0 (solid curve) and 
TV = 1.5 (dotted curve). 



FIG. 6. The behavior of the eigenfunctions of the m = 2 r-mode for the TV = 1 polytrope with r p /r e = 0.99 where r p andr e 
are the polar radius and the equatorial radius of the surface, respectively. Distributions along the redirection of the perturbed 
quantities are shown: 8u r (upper left panel), Svg (upper right panel), 5w v (lower right panel), and 8p (lower left panel). It is 
noted that the surface is located at r* = 1.0. In each panel, distributions of different 9 values are shown by different curves: A 
curve with No. L is the distribution at 9 = 7r/10 * L. 



FIG. 7. Same as Fig. 6 but for the rapidly rotating model with r p /r e = 0.71. The angular velocity is about 90% of that of 
the Kepler limit. 

FIG. 8. Critical rotational frequency for the m = 2 r-mode against the stellar temperature is plotted by the solid line for 
TV = 1 polytropes. This curve is obtained from the condition that l/r r = (see the text for details). This solid line divides the 
space into two distinct regions: one for the stable states and the other for the unstable states. In addition, the evolutionary 
curve is also plotted by the dashed line. This curve is obtained by equating the dissipation time scale and the cooling time 
scale, i.e r r = r coo i. The cooling time scale is assumed to be due to the modified URCA process. 



FIG. 9. Same as Fig. 11 but for the m = 3 mode. 



FIG. 10. Same as Fig. 11 but for the m = 4 mode. 



FIG. 11. Same as Fig. 11 but for the TV = 1.5 polytrope. 



FIG. 12. Same as Fig. 11 but for the TV = 0.5 polytrope. 

FIG. 13. Several evolutionary curves are shown for different cooling mechanisms. Meanings of the curves are shown in the 
figure. The thick solid curve is the critical curve of the TV = 1 polytrope for the m — 2 mode. It should be noted that the 
smallest values of the frequency of evolutionary curves are all located around 100Hz. 
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